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'How luscious lies the pea within the pod. . . " 

Emily Dickinson 



We evaluate the equation of state of the quasi-one-dimensional (ID) phase of C60 molecules in 
small carbon nanotubes, nicknamed "peas in a pod". The chemical potential and ID pressure are 
evaluated as functions of the temperature and density, initially with the approximation of near- 
est neighbor interactions and classical statistical mechanics. Quantum corrections and long-range 
interaction corrections are discussed, as are the effects of interactions with neighboring peapods. 
Transition phenomena involving the 3D coupling are evaluated. 



I. INTRODUCTION 



Carbon nanotubes have been found to exhibit a 
wide range of remarkable properties, concomitants of 
their nanometer cross-section and macroscopic length. 
Extensive research has recently explored the behav- 
ior of gases adsorbed within, outside of and in- 
between such nanotubesii 2 * 3 ^ 4 *^*^* 7 ". Among the most 
intriguing systems is that created when C@o bucky- 
balls lie within a tube of small radius R, since this 
is surely a quasi-one-dimensional (ID) system^. This 
system, nicknamed "peas in a pod" has begun to 
be studied using various theoretical and experimental 
methodaV ! 11 ! 12 ! 13 ! 14 ! 15 ! 16 ! 17 ! 18 . In this paper, we explore 
the equation of state of this system. Specifically, we eval- 
uate the ID pressure and the chemical potential, as func- 
tions of temperature and density. The calculation is rel- 
atively simple in the small tube case considered in this 
paper (e.g. (10,10) tubes, with R w 7 A) because the Ceo 
molecules are confined to the vicinity of the axis and the 
motion perpendicular to the axis can be factored out of 
the problem, as we shall see. The more general case in- 
volving tubes of varying size, with the attendant appear- 
ance of 3D effects, requires a more detailed study than 
is described here. Some numerical studies and molecular 
simulations exploring this problem have been performed 
by Hodak and Girifalco 14 i 15 i 16 . 



The next section discusses the model and presents the 
equation of state resulting from the simplest (classical, 
quasi-lD) procedure we have used. Section III assesses 
quantum and long-range corrections to this simple model 
and estimates the regime of applicability of those results. 
Section IV describes the effects of interactions with buck- 
yballs in neighboring tubes, including phase transitions. 



II. MODEL CALCULATIONS 

Girifalco and coworkers have showni&i^ how Cgo 
molecules within nanotubes of small radius (R f=a 7 A) 
experience potentials confining them to the vicinity of 
the tube axis. This paper employs the potentials devel- 
oped by that group to determine the equation of state of 
the quasi-lD buckyball fluid. We assume that the total 
potential energy Vtotai of the molecules satisfies 

V tota i{ Pi , Zi ) = J2u(pi) + V (\*i - z i\) (!) 

i i<j 

Here, = (p il Zi) is the position of the i-th molecule's 
center, written in terms of its coordinate zi parallel to the 
tube axis and its transverse displacement from the axis, 
P;. The potential U(p) is the molecule-tube interaction 
that determines the molecules' degree of localization near 
the axis, which depends also on the temperature T. The 
other interaction in Eq. (1), V(z), is that between a pair 
of Ceo molecules separated a distance z along the tube 
axis. Girifalco et al have shown that both interactions, 
U(p) and V(z), can be expressed in a universal way by 
the following reduced function^ 3 .: 

$(r) 

= U ^ V-gf ^ V° (2) 
3 V3.13f + 0.28/ 3 V3.13f + 0.28y w 

where ^(ro)! is the well depth at the equilibrium spacing 
ro, and f = (r — ri)/(ro — r\). The distance parameter 
r\ for each system (Ceo-Cgo, Ceo-tube) is chosen so that 
the interactions fit the universal function (2). Table I 
indicates the values of the used parameters. The func- 
tion in Eq. (2) is obtained using a continuum model for 
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the carbon surfaces. The carbon-carbon interatomic po- 
tential is integrated over the interacting surfaces assum- 
ing a mean surface density of carbon atoms and different 
Lennard-Jones parameters for atoms in Cgo-Cgo and C@o- 
graphene systems. Therefore, we ignore any dependence 
of the interactions on the atomic structure of either the 
tubes or the Ceo molecules. Similarly the neglect of any 
interactions that depend on the Cgo orientation means 
that the molecular rotational problem is not affected by 
the environment. This means that the rotational degree 
of freedom may be omitted completely from the present 
analysis 7 . 

Concerning the Ceo-Ceo interaction, there has been 
some discussion in the literature of the value of the 
asymptotic van der Waals coefficien1si2ii£ {Cq, the coeffi- 
cient of r~ e ), with values varying over the range 16.6 N 2 t 
eV A 6 to 29.05 N% t eV A 6 . The value used in Girifalco's 
calculations^ (Eq. (2)) is 20.0 N% t eV A 6 . 

In reality, the intermolecular interaction V ought to 
depend on the full 3D separation. By neglecting this in- 
teraction's dependence on p i and p^, Eq. (1) implicitly 
assumes that the molecules' rms transverse displacement 
d = \J (pf) is small compared to the mean separation 
a = L/N along the axis. We discuss the regime of valid- 
ity of that approximation below. Also implicit in Eq. (1) 
is the assumption that both the nanotube and the Ceo 
molecules are rigid, since otherwise their shape fluctua- 
tions would affect the nanotubes' potential energy. While 
conventional in this field, this assumption deserves fu- 
ture attention in the case of such strongly interacting 
molecules as Ceo- 

With Eq. (1), the properties of the system may be 
evaluated in terms of independent transverse and longi- 
tudinal contributions to the Hamiltonian. For example, 
the partition function Z of the Cgo molecules factorizes 
into transverse and longitudinal factors, Z = Z t Z Zl so 
that the Helmholtz free energy F is an additive function 
of these separate parts of the problem: F = F t + F z . We 
may therefore write the chemical potential p of the set 
of N molecules as a sum of a transverse part fi t and a 
longitudinal part /i z : 



U(p) = U + -p 2 + 7 P 4 



(5) 



The resulting transverse chemical potential is given 
by perturbation theory, assuming that the value of 7 is 
small: 



Pi 



P 



— coth (-} 
mui V 2/ 



(6) 



Here x = (3Ttio with to = yjk/m, and the regime of 
validity of this expression is that range of T for which 
the rightmost term in fi t is small compared to the middle 
term (i.e., k B T < k B T± = k 2 /j). In the case of (10,10) 
nanotubes, the upper limit to the T range for which this 
expression is valid is about 13,000 K (based on our fits 
to the Girifalco potential U(p): k = 1.32 eV/ A 2 and 
7 = 1.58 A~ 4 ). Thus, the expansion suffices in practice 
throughout the relevant range of T. 

The present case corresponds to x « 32 K/T, which 
is much less than 1 at room temperature. Therefore, it 
is appropriate to substitute the fully classical expression 
for n t : 



Ppt 



1» 1 W 

Arp 



dh 



-fiU(p) 



(7) 



Here At = {2tt(3Ti 2 /m) 1 / 2 is the de Broglie thermal 
wavelength of the molecule, which is about 0.04 A at 
room temperature. The resulting expression for fi t in the 
case of a harmonic potential with a small quartic term 
coincides with the classical limit (a; -C 1) of the preced- 
ing quantum expressions for both \i t and the transverse 
contribution to the specific heat (per molecule), C t /N: 



Pt = U + - ln(ir) + 



8-, 



{(3ui 2 m) 2 



(8) 



/i = fi t + p z 



(3) 



C t 
Nk f 



= 2 ^ = 2-16^ 

/?(w 2 m) 2 T 4 



(9) 



One may immediately evaluate pt from the transverse 
partition function Z t because that part of the problem 
consists of independent motions of individual molecules: 



&IH = -jy In Z t = In 



-/3E a 



(4) 



Here = k B T and the sum is over all eigenvalues E a 
of the Schrodinger equation for transverse motion of a 
single molecule. This equation may be evaluated analyt- 
ically in the case when the transverse potential is nearly 
harmonic, i.e. when U includes just a small quartic term, 
involving a parameter 7: 



This includes a constant contribution (expected for two 
harmonic motions in the classical regime of T), plus an 
anharmonic term proportional to T . 

Figure 1 shows p t as a function of T, computed using 
the different aproximations described above. The full line 
is the quantum result (Eq. (6)) obtained assuming the 
nearly harmonic transverse potential of Eq. (5). It is 
interesting to note that the corresponding classical ex- 
pression, Eq. (8), is extremely similar to the quantum 
one (they coincide in the figure), showing that quantum 
effects are really negligible. The dashed line shows the 
classical result Eq. (7) using the exact potential. For 
T > 500 K, the curves start to differ as the real shape of 
the potential deviates from the quartic aproximation. 
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The solution of the remaining ID problem of interact- 
ing buckyballs requires a numerical calculation, in gen- 
eral. Quantum calculations have been carried out for 
some ID fluids (helium and hydrogen 2 ^ in their ground 
states), but for classical systems at finite T the deter- 
mination of the equation of state is more straightfor- 
ward. Specifically, the equation of a ID classical system 
assumes an analytic form in the case when nearest neigh- 
bor interactions adequately represent the total potential 
energy 2 ^. We assume that to be the case here (and jus- 
tify that in Section III) so that we may write down an 
expression relating the ID density, N/L = 1/a, in terms 
of the ID pressure P: 

L = J Q do) 

a f °° dz z e-fivl*)+*P\ [ ' 

where V(z) is the Cg -C 6 o interaction given by Eq. (2). 
Similar results were reported recently by Hodak and 
Girifalcoi&it 

Figure 2 shows the "compression ratio" , i.e. the ratio 
of P to the pressure of an ideal ID gas (fc^T/a); the in- 
teraction effects can be clearly observed. The pressure is 
evidently reduced by interactions over an extended den- 
sity range until the point when the molecules start to 
come nearly into contact (a ~ 12.5 A at T — 300 K), 
above which density the pressure rises extremely rapidly. 
One observes that the Ceo gas behavior is nearly ideal at 
T > 1000 K while the effects of the interaction are quite 
significant below room temperature. The Boyle tempera- 
ture is around 1400 K. The calculated behavior and high 
characteristic temperature are not surprising in view of 
the fact that the well depth of the pair potential is about 
3200 K. 

The ID pressure can be integrated (according to the 
Gibbs-Duhem relation) to obtain the total chemical po- 
tential of the Ceo fluid at temperature T and density 1/a: 

fi{T,a) = fi t (T)+n id (T,a ) + [ a(P' 7 T)dP' (11) 

J Pa 

Here Hid(T, do) = fcsTln(AT/ao) is the chemical po- 
tential of a ID ideal gas at a reference density (1/oq) that 
is low enough for the fluid to be ideal^i. The integration 
domain extends from the reference pressure Pq = ksT/ao 
to the pressure corresponding to the specified spacing a. 
The result of this analysis appears in Figure 3 for a range 
of temperatures. 

Figure 4 exhibits the effects of the interaction on the 
chemical potential, expressed in terms of the difference 
between the chemical potential of the (fully interacting) 
system and that of an ideal system which has the mutual 
interactions "turned off" : 

A fii nt (T, a) = fi(T,a)~[fi t {T)+fi id (T,a)} = J (a(P',T) 

(12) 



The rightmost integral may alternatively be expressed 
in terms of the ID compressibility, k — — (d In a/dP)T as 
follows: 

Afi lnt (T,a) = f -7^) da ' ( 13 ) 

Joo V a K[a ,T)J 

Because the attraction dominates the effect of the in- 
teraction over most of the density range shown in Figures 
2 to 4, the chemical potential shift due to the interaction 
is negative until one reaches high density. 

III. CORRECTIONS TO THE MODEL 

We assess three corrections to the preceding section's 
results. One is the long range correction to the thermo- 
dynamic properties. This is the energy shift associated 
with the neglect of interactions other than nearest neigh- 
bors and could be written in terms of the probability per 
unit length P>(x) that a second, third or other (than 
first) nearest neighbor is present at distance x. Since we 
seek only an estimate of its effect, we take P>(x) to be a 
Heaviside step function for \x\ > 2a, where a is the hard 
core diameter (7.1 A) of the Ceo molecule. Then, our 
approximate energy correction per molecule due to this 
neglect is given by 

r°° dz 

AE > /N= / — V(z) (14) 

Note that a "missing" factor of (1/2) is cancelled by 
the omission of the region z < 0. If the potential V(z) 
were the familiar Lennard-Jones 6-12 interaction, the re- 
sult of this integration would be AE > /(Ne) = — a/ (40a) 
(where e is the well depth), which represents a negligible 
correction of the energy. 

The second correction we assess is the quantum cor- 
rection. To do this, we introduce a phonon model as a 
means of characterizing the dynamics of the fluid. As 
was shown recently^, for a Lennard-Jones system in its 
ground state, the ratio of the zero-point energy E zp to 
the pair potential well depth e is 

ElL = 2 4 / 3 -1 A* (15) 
e it 

Here A* = /i/fo^me) 1 / 2 ] is the de Boer quantum pa- 
rameter. In the case of a Cqq modeled by such an 
interaction 22 , A* w 0.003, so this ratio is 0.002. This 
indicates that the quantum effects are small even in 
the ground state. They become even smaller (relatively 
speaking) with increasing T and can be neglected at 
T > E Z p/kB, which is of order 10 K for Cgo- 
^Wm kll^pwe estimate the effect on the mutual interac- 
tioH ejlergy of the molecules' motion perpendicular to the 
axis. To do so, we evaluate the expectation value of the 
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difference in potential energy due to changing the approx- 
imate separation of a pair of molecules (zy ) used in Eq. 
(1) into their fully 3D separation r tj = [p 2 - + zfj) 1 / 2 - To 
lowest order in the transverse separation, the difference 
in potential energy is AV(z) = V'(z)[p 2 j/(2zij)], where 
the prime means a derivative. Substituting the mean 
value (/£>?•) = 2d 2 , the energy shift per molecule due to 
this substitution is obtained from the radial distribution 
function gio{z)'- 

= / dz g 1D {z) (16) 

For this estimate, we substitute for gn)(z) a step func- 
tion at z = a. The result for a Lennard- Jones potential 
is 



AF inter = / dz V[(b 2 + Z 2 ) 1 ' 2 ] (18) 

Here b is the separation between tubes and v is the 
number of nearest neighbor tubes. This expression as- 
sumes a simple form in the case when b is sufficiently 
large that the interaction may be approximated by the 
asymptotic form V(r) sa —Cs/r 6 , in which case the inte- 
gration yields 

AF mter = -vN (19) 
16 o° a 

Including this term in the equation of state yields a 
relation for the shift of the ID pressure proportional to 
1/a 2 



AE 1 = _2A^ 

Ne 91 aa K ' 

For a (10,10) tube in a high density case, with a = a ks 
9.5 A, for example, at room temperature (d = 0.5 A), 
we obtain AE t /(Ne) — 6xl0~ 4 , which is negligibly small 
(using e w 3200 K). Although this is a relatively crude 
estimate, the order of magnitude is plausible in view of 
the small ratio of d/a. Evidently, the small value of AE t 
implies that the ID approximation is adequate for the 
present circumstances. 



IV. TRANSITIONS DUE TO INTERACTIONS 
WITH NEIGHBORING TUBES 

In this section we evaluate an interesting possibility- 
transitions of the Ceo fluid as a result of molecular inter- 
actions between buckyballs in neighboring tubes. In our 
analysis, we assume an infinite array of perfectly parallel 
and identical tubes and chemical equilibrium among the 
Ceo molecules. Under these circumstances, molecules in 
nearby tubes can undergo a cooperative transition, a con- 
densation to a 3D liquid. A qualitatively similar problem 
(condensation of helium confined within interstitial chan- 
nels of a carbon nanotube bundle) was studied previously 
by our group, using a Monte Carlo evaluation of a lattice 
gas modei&. 

The present analysis is carried out first by perturba- 
tion theory; then, a comparison is made with an ex- 
actly solved model that captures the essential physics 
but employs a number of simplifications. Although nei- 
ther model is completely satisfactory, we will see that the 
two sets of results are quite similar because the key pre- 
dictions (parameters of the critical point) turn out to be 
only weakly dependent on the assumptions. 

The perturbation theory is based on the shift in free en- 
ergy of the system due to interactions between molecules 
in neighboring tubes, AFi nter : 



Note that a = (ivn / lQ)bV {b); thus it is proportional 
to the product of the strength of the transverse inter- 
action and the number of neighbors interacting with a 
molecule (~ b/a). 

Figure 5 shows the effect on the isotherms of including 
the intertube interaction term. One observes that van 
der Waals loops develop, indicating the appearance of a 
phase transition. Figure 6 shows how the critical temper- 
ature T c varies with the parameter a. For a typical lattice 
of (10,10) tubes with b = 17 A and v = 6, a « 2180 K 
Awhich yields T c w 530 K. There is a singular rise of 
T c for small a, followed by a region of weaker depen- 
dence. We may understand this behavior by examining 
an alternative model, the anisotropic Ising model applied 
to the lattice gas version of this problem. As shown by 
Fisher—, for the case of such a model with nearest neigh- 
bor interactions J/ along the ID chains and transverse 
interactions Jt, the critical temperature is given by the 
expression 



B c ~ ln(l/c) - HHl/c)} [ ' 

where c = Jt/Ji- This is an asymptotic formula valid 
when c < 0.1 (weak transverse interaction). Figure 6 
shows the dependence derived from this expression when 
the values <// = 800 K and J t = 18.2 K are substituted 
in the Fisher formula. The result is seen to track the 
perturbation theory behavior well. This similarity is sur- 
prising, at first sight, because the perturbation theory 
is a mean field theory, neglecting density fluctuations in 
adjacent tubes. We rationalize this agreement by noting 
that the exact calculation yields critical exponents (3D 
Ising model) that are, in fact, different from those of the 
mean field perturbation theory. On the other hand, in 
3D, one does not expect a large error in the mean field 
theory value of the critical temperature^. 
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System 


|$(r )| (eV) 


ro (A) 


n (A) 


C6o-Ceo 


0.278 


10.05 


7.10 


C 60 -(10,10) tube 


0.537 


13.28 


10.12 



TABLE I: Parameters of the universal function Eq. (2) for 
the carbon structures investigated. 



where S = b/a. The difference would be identically zero 
if the Euler-MacLaurin theorem were analitically used to 
evaluate it. The integral can be analitically evaluated 
yielding 37r/(8<5 5 ). Hence, we compute the sum numeri- 
cally and find the difference. It is negligible small when 
b/a > 2 and increases as b/a goes to zero. This is shown 
in Figure 7, expressed as a freezing temperature 



Finally, we note that Carraro has described the possi- 
bility of a crystallization transition for such a problem 24 . 
To estimate the onset temperature for this freezing, we 
simply evaluate the energy difference AEf between the 
ground state crystal and the energy of the alternative 
fluid state. This is determined by comparing the energy 
computed above (AFi nter ) for the intertube interaction 
energy. The difference is just the difference between the 
integral and the sum that arises when the buckyballs are 
located in lattice sites, separated by a: 



AE f = - 



vC 6 
2a 6 



E 



(<5 2 + n 2 ) 3 



dz 



{S 2 + z 2 f 



(22) 



T f = AE f /k B 



(23) 



In the case of Ceo inside (10,10) tubes, b/a = 1.7 and 
Tf ~ 0.5K. Such small values of Tf imply that such a 
transition will be difficult to observe in the laboratory. 
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FIG. 1: The transverse chemical potential fit as function of temperature. The full line corresponds to the quartic aproximation 
of the potential (Eq.5); quantum and classical results (Eqs. (6) and (8)) are so similar that they cannot be distinguished over 
the whole range of T. The dashed line shows the classical result using the full potential, Eq. (7). 




FIG. 2: Compression ratio at various temperatures of the ID pressure P to that of an ideal ID gas, Pideai = 
function of the density 1 /a. 
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FIG. 3: Total chemical potential /i(T, a) = nt(T) + /i z (T, a) of the C60 fluid as a function of molecular density 1/a and 
temperature T. 




FIG. 4: Chemical potential shift due to the intermolecular interactions A/zi n t(T, a) as a function of 1/a. 




FIG. 5: Effects of the interaction between C60 moiecuies in nearest neighbor tubes (a — 2180 K A). The isotherms develop van 
der Waals loops as a function of the temperature, with a critical temperature around 530 K. 




FIG. 6: Dependence of the critical temperature T c on the strength of the transverse interaction and geometry of the tube 
lattice, given by the parameter a. The full curve corresponds to the perturbation theory results and the dashed one indicates 
the anisotropic lattice model predictions (see text). 
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FIG. 7: Energy difference AEf between the ground state crystal and the fluid state as a function of the dimensionless separation 
between tubes 8 — b/a. 



